function [pX,pY,P2d,ODraw] = ProcessImageGUI_v2(handles)

% [pX,pY,OD] = ProcessImage(hdfFilePath,autoROI,userROI,doFitData)
% Finds Gaussian-shaped signals on fluorescence or absorption images.
% Returns:
% pX - The parameters for a Gaussian fit to the data in the x (vertical)
% direction.
% pX(1) = Gaussian amplitude
% pX(2) = Peak position (in units of image pixels)
% pX(3) = Standard deviation (in units of image pixels)
% pX(4) = Offset
% pY - Same description as pX
% OD - The image (OD if it's an absorption image, raw image for a fluor.)
% num - Counts in raw "atoms" image, with electronic noise subtracted
% den - Counts in raw "background" image, with electronic noise subtracted.

fit2d=handles.do2DFit;

% Inputs from hdf5 file
rawimages = cast(h5read(hdfFilePath,'/Image/image'),'double');
img_var = h5read(hdfFilePath,'/Image/configuration');
img_time = h5read(hdfFilePath,'/Image/exposuretime');
type_img = h5read(hdfFilePath,'/Image/typeimaging');
CCDbin = 1; % Amount of binning at camera sensor; '1' is 1x1 binning and '2' is 2x2 binning

%%% TODO: Check validity of image, has pixels, correct size, has been
%%% binned etc.

%Software binning
bindim=[1,1];  % Software binning size
rawimages=BinImg2d(rawimages,bindim);

% Libraries used in this function
expp = ParamImgSys(img_var);
cons = Constants();

dim = size(rawimages)  % size of hdf5 image group
rows = dim(1);          % rows in one image file (# of horiz. pixels)
cols = dim(2);          % columns in one image file (# of vert. pixels)


ovflw = 'no';   % setting 'ovflw' flag to 'no'
if  type_img == 3 % Absorption  %%size(dim,2)
    elec = rawimages(:,:,3);
    num = rawimages(:,:,1) - elec;
    den = rawimages(:,:,2) - elec;
    maxcntN = max(max(num));
    maxcntD = max(max(den));
    if maxcntN > 2^16 || maxcntD > 2^16
        ovflw = 'yes';
        warning('Bin overflow warning');
    end
    
    %scaling the baselevel of the atom image
    numforscale=10;  %using a 10x10 pixel area to determine scaling fraction
    numleftup = sum(sum(num(1:numforscale,1:numforscale)));
    numleftdn = sum(sum(num(rows-numforscale:rows,1:numforscale)));
    numrightup = sum(sum(num(1:numforscale,cols-numforscale:cols)));
    numrightdn = sum(sum(num(rows-numforscale:rows,cols-numforscale:cols)));
    baselevelnum = numleftup + numleftdn + numrightup + numrightdn;
    
    denleftup = sum(sum(den(1:numforscale,1:numforscale)));
    denleftdn = sum(sum(den(rows-numforscale:rows,1:numforscale)));
    denrightup = sum(sum(den(1:numforscale,cols-numforscale:cols)));
    denrightdn = sum(sum(den(rows-numforscale:rows,cols-numforscale:cols)));
    baselevelden = denleftup + denleftdn + denrightup + denrightdn;
    blscale = (baselevelden/baselevelnum); %resulting baselevel scaling
    modnum=num*blscale;
    ODraw = -log((num./den)');  %SIMON ROTATES THE MATRIX HERE: (num./den)'
    ODrawmod = -log(modnum./den);
    
    Alow = 0;   % Amplitude lower limit for fit
    Ahi = 5;    % Amplitude upper limit for fit
elseif type_img == 1 % Fluorescence
    ODraw = rawimages';  %SIMON ROTATES THE MATRIX HERE
    ODraw = ODraw - min(min(ODraw));
    
    Alow = 10;   % Amplitude lower limit for fit
    Ahi = 65000;    % Amplitude upper limit for fit
else
    error('Can''t determine imaging method, type_img error.');
end

% Filter image to find ROI
N = 60;  %amplitude of 2D gaussian
sig = 30;  %size of 2D gaussian
ODfilt = conv2(ODraw,gaussian2d(N,sig),'same'); % Low-pass filter image

yC = 1:expp.short_px_nb;   %y-axis IS VERT. DIMENSION OF SENSOR/HORZ. dim. of expt.
xC = 1:expp.long_px_nb;

autoROIok = false;

if autoROI % Find ROI and extract first approx. signal from filtered image
    % Determine pre-ROI for fitting (for initial guesses)
    mx = max(max(ODfilt));
    [I,J] = find(ODfilt == mx);
    wndw = 50;
    if ~(I - wndw < 1 || I + wndw > max(xC) || J - wndw < 1 || J + wndw > max(yC))
        autoROIok = true;
        
        ROI(1) = I(1);
        ROI(2) = J(1);
        
        % Projections of ODfilt using pre-ROI as size (for initial guesses)
        yIntFilt = sum(ODfilt(:,ROI(2)-wndw:ROI(2)+wndw),2);
        xIntFilt = sum(ODfilt(ROI(1)-wndw:ROI(1)+wndw,:),1);
        
        % Least-squares fit of ODfilt (initial guesses)
        ampGuessY = max(yIntFilt) - min(yIntFilt);
        ampGuessX = max(xIntFilt) - min(xIntFilt);
        options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
        pYtmp = fminsearch(@(p)FitGauss(p,yC,yIntFilt),[ampGuessY ROI(1) 100 min(yIntFilt)],options);
        pXtmp = fminsearch(@(p)FitGauss(p,xC,xIntFilt'),[ampGuessX ROI(2) 100 min(xIntFilt)],options);
        
        % Determining autoROI based on size of ODfilt (for real fit)
        sigmas = 3;
        xWndw = round(sigmas*pXtmp(3));
        yWndw = round(sigmas*pYtmp(3));
        xROI = ROI(2)-xWndw:ROI(2)+xWndw;
        yROI = ROI(1)-yWndw:ROI(1)+yWndw;
        
        % Remove parts of the ROI if it goes outside of the photo.
        if min(xROI) < 1
            xROI = xROI(xROI > 0);
        end
        if min(yROI) < 1
            yROI = yROI(yROI > 0);
        end
        
        if max(yROI) > expp.short_px_nb
            yROI = yROI(yROI < expp.short_px_nb + 1);
        end
        if max(xROI) > expp.long_px_nb
            xROI = xROI(xROI < expp.long_px_nb + 1);
        end
        
        % Projections of ODraw using autoROI as size (for real fit)
        yInt = sum(ODraw(yROI,xROI),2);
        xInt = sum(ODraw(yROI,xROI),1);
        
    else
        warndlg('Can''t auto-choose ROI: Estimated signal close to edge of image')
        autoROIok = false;
    end
end

if ~autoROIok
    % Use user-input ROI parameters to guess initial fitting parameters
    % User-defined ROI
    yROI = userROI(3):userROI(4);
    xROI = userROI(1):userROI(2);
    
    % Projections of ODraw using user-defined ROI as size (for real fit)
    yInt = sum(ODraw(yROI,xROI),2);
    xInt = sum(ODraw(yROI,xROI),1);
    
    % Amp. and sizes (initial guesses)
    pYtmp(1) = max(yInt) - min(yInt);
    pXtmp(1) = max(xInt) - min(xInt);
    pYtmp(2:3) = [mean(userROI(3:4)) 150];
    pXtmp(2:3) = [mean(userROI(1:2)) 150];
end

if doFitData % Real fit of ODraw using initial guesses
    options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
    pY = fminsearch(@(p)FitGauss(p,yROI,yInt),[pYtmp(1:3) min(yInt)],options);
    pX = fminsearch(@(p)FitGauss(p,xROI,xInt'),[pXtmp(1:3) min(xInt)],options);
else
    pY = zeros(1,4);
    pX = pY;
end


if fit2d==1  %2d fitting procedure
    % Generating guesses and bounds by smoothing OD
    smoothsize=8;
    filtermat=ones(smoothsize,smoothsize);
    ODsmooth=filter2(filtermat,ODraw(:,:))/smoothsize^2; %ODraw, smoothed alot for finding guesses
    
    figure(100)
    subplot(2,3,1)
    surf(ODsmooth,'Marker','.') %Plot smoothed data
    cax=caxis;
    colorbar
    shading flat
    grid off
    title('Smoothed ODraw')
    view(0,90);
    xlabel('x')
    ylabel('y')
    
    amp2d=max(max(ODsmooth));
    [yoff2d,xoff2d] = find(ODsmooth == amp2d);
    sigx2d=rows/10;
    sigy2d=cols/10;
    guessave=5; %using a 5x5 pixel area to determine scaling fraction
    guessleftup = sum(sum(ODsmooth(1:guessave,1:guessave)))/guessave^2;
    guessleftdn = sum(sum(ODsmooth(cols-guessave:cols,1:guessave)))/guessave^2;
    guessrightup = sum(sum(ODsmooth(1:guessave,rows-guessave:rows)))/guessave^2;
    guessrightdn = sum(sum(ODsmooth(cols-guessave:cols,rows-guessave:rows)))/guessave^2;
    offset2d=(guessleftup + guessleftdn + guessrightup + guessrightdn)/4 ;
    slopex2d=((guessrightup + guessrightdn) - (guessleftup + guessleftdn))/(2*cols);
    slopey2d=((guessleftdn + guessrightdn) - (guessrightup + guessleftup))/(2*rows);
    
    IniCondNR=[offset2d,amp2d,xoff2d,yoff2d,sigx2d,sigy2d,slopex2d,slopey2d];
    
    % Creating matrices for fitting
    xvec=1:rows;
    yvec=1:cols;
    [xdata,ydata]=meshgrid(xvec,yvec);
    ODerror2d(1:rows*cols)=1;    % DUMMY
    numberofpoints=cols*rows*ones(rows,cols); % matrix in which each entry is the number of points
    fitcoordmatrix=[xdata(:),ydata(:),ODerror2d(:),numberofpoints(:)];
    
    Z2d=ODraw(:);  %flattened(1D) vector of unsmoothed image
    % Fitting - data is treated with equal weight
    [P2d,r2d,J2d]=nlinfit(fitcoordmatrix,Z2d(:),@gauss2derror,IniCondNR);
    
    % Generate OD plot based on fit parameters
    dummyerrorod(1:rows*cols)=1;
    dummynumberofpoints=ones(rows,cols);
    dummycoordmatrix=[xdata(:),ydata(:),dummyerrorod(:),dummynumberofpoints(:)];
    ODcalcvec=gauss2derror(P2d,dummycoordmatrix);           %OD calculated from fit paramters
    ODcalc = reshape(ODcalcvec, size(xdata));
    ODcalcICvec=gauss2derror(IniCondNR,dummycoordmatrix);   %OD based on initial conditions
    ODcalcIC = reshape(ODcalcICvec, size(xdata));
    ODresiduals = ODsmooth-ODcalc;
    
    subplot(2,3,3)
    surf(ODcalcIC,'Marker','.')  %Plot 2d fit using IC
    %caxis(cax)
    colorbar
    shading flat
    grid off
    view(0,90);
    xlabel('x')
    ylabel('y')
    title('OD from initial conditions')
    hold on
    
    subplot(2,3,2)
    surf(ODcalc,'Marker','.')  %Plot 2d fit
    caxis(cax)
    colorbar
    shading flat
    grid off
    view(0,90);
    xlabel('x')
    ylabel('y')
    title('OD from 2d fit')
    hold on
    
    %     % Removing background and slope from raw data
    %     bgslopefit=(xdata(:)-P2d(3))*P2d(7)+(ydata(:)-P2d(4))*P2d(8)+P2d(1);
    %     ODcorrect=ODraw(:)-bgslopefit(:);
    %     ODcorrect2d=reshape(ODcorrect,size(xdata));
    
    subplot(2,3,4)
    surf(ODresiduals,'Marker','.')  %Plot 2d fit
    caxis('auto')
    view(0,0);
    colorbar
    shading flat
    grid off
    title('True residues (NOT weighted)')
    xlabel('x')
    ylabel('y')
    hold off
    
    subplot(2,3,5)
    surf(ODresiduals,'Marker','.')  %Plot 2d fit
    caxis('auto')
    %view(0,0);
    colorbar
    shading flat
    grid off
    title('True residues (NOT weighted)')
    xlabel('x')
    ylabel('y')
    hold off
end %2d fitting procedure

% Give the width as the std deviation
widthYm = (P2d(6))*expp.px_size*CCDbin/expp.mag;    %Y cloud size (m)
widthXm = (P2d(5))*expp.px_size*CCDbin/expp.mag;    %X cloud size (m)
amp = (pX(1) + pY(1))/2                    % "2D" Gaussian fit amplitude (mean of single Gaussian fits)

if type_img == 1 % Fluorescence
    %     % Estimate number of atoms using camera
    %     capture_fraction = 0.5*(1 - sqrt(1 - (expp.na_telec)^2)); % fraction of total light captured by lens ( ~ 0.5 %)
    %     total_num_electrons = amp*sqrt(2*pi)*mean([pX(3) pY(3)]);
    %     total_photons = total_num_electrons/(expp.qe_std*expp.trans780*expp.trans_telec*capture_fraction);
    %     num_atoms = 2*total_photons/(img_time*cons.Gamma) % num_atoms = (num_photons/time) / (scattering rate)
    
    % Estimate number of atoms using camera
    total_num_electrons = P2d(2)*sqrt(2*pi)*mean([P2d(6) P2d(5)]);
    total_photons = total_num_electrons/(expp.qe*expp.coll_eff);
    num_atoms = 2*total_photons/(img_time*1e-6*cons.Gamma) % num_atoms = (num_photons/time) / (scattering rate)
    set(handles.text_num_atoms,'string',['# Atoms: ' num2str(num_atoms,'%.2e')]);
elseif size(dim,2) == 3 % Absorption
    totalODraw = sum(sum(ODraw))*(expp.px_size*CCDbin)^2;
    num_atomsOD = totalODraw/AbsCross(0,0);
    num_atomsFit = P2d(2)*2*pi*widthYm*widthXm/AbsCross(0,0);
    set(handles.text_num_atoms,'string',['# Atoms (img): ' num2str(num_atomsOD,'%.2e')]);
end

% if type_img == 1 % Fluorescence
%     figure(1)
%     subplot(3,3,[1 2 4 5])
%     oneMat = ones(size(ODraw));
%     zeroMat = nan(size(ODraw));
%     zeroMat(yROI,xROI) = 1;
%     ODwindow = ODraw.*zeroMat;
%     imagesc(ODwindow);
% imagesc((ODfilt));
axes(handles.axes_processed)
imagesc(ODraw)
line([xROI(1) xROI(end)],[yROI(1) yROI(1)])
line([xROI(1) xROI(1)],[yROI(1) yROI(end)])
line([xROI(end) xROI(end)],[yROI(1) yROI(end)])
line([xROI(1) xROI(end)],[yROI(end) yROI(end)])
axis equal
xlims = xlim;
ylims = ylim;
set(gca,'FontSize',10)

% subplot(3,3,[3 6])
axes(handles.axes_y_proj)
plot(yInt,yROI,gauss(pY,yC),yC);
set(gca,'YDir','reverse');
ylim(ylims);
set(gca,'FontSize',10)

% subplot(3,3,[7 8])
axes(handles.axes_x_proj)
plot(xROI,xInt,xC,gauss(pX,xC));
xlim(xlims);
set(gca,'FontSize',10)

if type_img == 3 % Abs. img.
    axes(handles.axes_raw_abs_atoms)
    imagesc(y(:,:,1))
    axis off
    
    axes(handles.axes_raw_abs_no_atoms)
    imagesc(y(:,:,2))
    axis off
    
    axes(handles.axes_raw_background)
    imagesc(y(:,:,3))
    axis off
end

set(handles.text_x_amplitude,'string',['X amp.: ' num2str(pX(1),'%.2e')]);
set(handles.text_y_amplitude,'string',['Y amp.: ' num2str(pY(1),'%.2e')]);
set(handles.text_x_width,'string',['X width: ' num2str(pX(3),'%.2e')]);
set(handles.text_y_width,'string',['Y width: ' num2str(pY(3),'%.2e')]);

%     h = subplot(3,3,9,'Visible', 'off');
%     cla(h);
%     text(0, .5, sprintf('%s\n%s\n%s\n%s', ...
%         ['v-width (um): ' num2str(widthXm*cons.um)], ['h-width (um): ' num2str(widthYm*cons.um)], ...
%         ['peak(v,h) = (' num2str(round(pX(2))) ',' num2str(round(pY(2))) ')'], ...
%         ['Na = ' num2str(num_atoms,'%e')]), ...
%         'VerticalAlignment', 'middle', ...
%         'FontName', 'Courier New', ...
%         'FontWeight', 'bold', ...
%         'FontSize', 12);

%     elseif type_img == 3 % Absorption
%     figure(1)
%     subplot(4,3,[1 2 4 5 ])
%     oneMat = ones(size(ODraw));
%     zeroMat = nan(size(ODraw));
%     zeroMat(yROI,xROI) = 1;
%     ODwindow = ODraw.*zeroMat;
%     imagesc(ODwindow);
%     % imagesc((ODfilt));
%     line([xROI(1) xROI(end)],[yROI(1) yROI(1)])
%     line([xROI(1) xROI(1)],[yROI(1) yROI(end)])
%     line([xROI(end) xROI(end)],[yROI(1) yROI(end)])
%     line([xROI(1) xROI(end)],[yROI(end) yROI(end)])
%     axis equal
%     xlims = xlim;
%     ylims = ylim;
%
%     subplot(4,3,[3 6])
%     plot(yInt,yROI,gauss(pY,yC),yC);
%     set(gca,'YDir','reverse');
%     ylim(ylims);
%
%     subplot(4,3,[7 8])
%     plot(xROI,xInt,xC,gauss(pX,xC));
%     xlim(xlims);
%
%     subplot(4,3,10)
%     imagesc(num')
%     axis equal
%
%     subplot(4,3,11)
%     imagesc(den')
%     axis equal
%
%     subplot(4,3,12)
%     imagesc(elec')
%     axis equal
%
%     h = subplot(4,3,9,'Visible', 'off');
%     cla(h);
%     text(0, .6, sprintf('%s\n%s\n%s\n%s\n%s\n%s \n%s \n%s \n%s', ...
%         ['v-width (um): ' num2str(widthXm*cons.um)], ['h-width (um): ' num2str(widthYm*cons.um)], ...
%         ['peak(v,h) = (' num2str(round(pX(2))) ',' num2str(round(pY(2))) ')'], ...
%         ['N(image) = ' num2str(num_atomsOD,'%e')], ...
%         ['N(fit) = ' num2str(num_atomsFit,'%e')], ...
%         ['peak pixel cnt = ' num2str(max(max(den)))], ...
%         ['Peak OD(image) = ' num2str(max(max(ODraw)))], ...
%         ['Peak OD(fit) = ' num2str(P2d(2))], ...
%         ['Mean OD(image) = ' num2str(mean(mean(ODraw)))]), ...
%         'VerticalAlignment', 'middle', ...
%         'FontName', 'Courier New', ...
%         'FontWeight', 'bold', ...
%         'FontSize', 12);
% end

% figure(3)
% imagesc(rawimages(:,:,3))
% axis equal

% cnt = numel(rawimages(:,:,2));
%
% sigAt = reshape(rawimages(:,:,1),1,cnt); % signal with atoms
% sigBK = reshape(rawimages(:,:,2),1,cnt); % background signal (no atoms)
%
% [Nat,pAt] = hist(sigAt,100);
% [Nbk,pBK] = hist(sigBK,100);
%
% figure(2)
% plot(pAt,Nat,pBK,Nbk)
%
% figure(3)
% imagesc(rawimages(:,:,1))